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Abstract: Complex colloidal fluids, such as emulsions stabilized by particles with 
complex shapes, play an important role in many industrial applications. However, 
understanding their physics requires a study at sufficiently large length scales while 
still resolving the microscopic structure of a lai^ge number of particles and of the lo- 
cal hydrodynamics. Due to its high degree of locality, the lattice Boltzmann method, 
when combined with a molecular dynamics solver and parallelized on modern super- 
computers, provides a tool that allows such studies. Still, mnning simulations on hun- 
dreds of thousands of cores is not trivial. We report on our practical experiences when 
employing large fractions of an IBM Blue Gene/P system for our simulations. Then, 
we extend our model for spherical particles in multicomponent flows to anisotropic 
ellipsoidal objects rendering the shape of e.g. clay particles. The model is applied to 
a number of test cases including the adsorption of single particles at fluid interfaces 
and the formation and stabilization of Pickering emulsions or bijels. 

Keywords: Complex Colloidal Fluid, Lattice Boltzmann method. Blue Gene/P, 
Domain Decomposition, Parallelization, Ellipsoid. 



1. Introduction 

Colloidal particles are highly attractive in the food, cosmetics, and medical industries to stabilize 
emulsions or to develop sophisticated ways to deliver drugs at the right position in the human body. The 
underlying microscopic processes of emulsion stabilization with particles can be explained by assuming 
an oil-water mixture. Without additives, both liquids phase separate, but the mixture can be stabilized 
by adding small particles which diffuse to the interface and stabilize it due to a reduced interfacial 
free energy. If for example individual droplets of one phase are covered by particles, such systems ai^e 
referred to as "Pickering emulsions", which have been known since the beginning of the 20th century |[l|, 
0]. Particularly interesting properties of such emulsions are the blocking of Ostwald ripening and the 
rheological properties due to irreversible particle adsorption at interfaces or interface bridging due to 
particle monolayers liSJ. Recently, interest in paiticle-stabilized emulsions has led to the discovery of 
a new material type, the "bicontinuous interfacially jammed emulsion gel" (bijel), which shows an 
interface between two continuous fluid phases that is covered by particles. The existence of the biiel 
was predicted in 2005 by Stratford et al. p}] and experimentally confirmed by Herzig et al. in 2007 Js]. 

Computer simulations are promising to understand the dynamic properties of particle-stabilized mul- 
tiphase flows. However, the shortcomings of traditional simulation methods quickly become obvious: a 
suitable simulation algorithm is not only required to deal with simple fluid dynamics but has to be able to 
simulate several fluid species while also considering the motion of the particles and the fluid-particle in- 
teractions. Some recent approaches trying to solve these problems utilize the lattice Boltzmann method 
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for the description of the solvents i^. The lattice Boltzmann method can be seen as an alternative to 
conventional Navier-Stokes solvers and is well-established in the literature. It is attractive for the cur- 
rent application since a number of multiphase and multicomponent models exist which are comparably 
straightforward to implement. In addition, boundary conditions have been developed to simulate sus- 
pended finite-size particles in flow. These ai^e commonly used to study the behavior of particle-laden 
single phase flows A few groups combined multiphase lattice Boltzmann solvers with the known 
algorithms for suspended particles |l,[8l]. In this paper we follow an alternative approach based on the 
multicomponent lattice Boltzmann model of Shan and Chen [QJ which allows the simulation of multiple 
fluid components with surface tension. Our model generally allows arbitrary movements and rotations 
of rigid particles of arbitrary shape. Further, it allows an arbitrary choice of the particle wettability - 
one of the most important parameters for the dynamics of multiphase suspensions [3]. For a detailed 
introduction to the method see Ref. M\, where our model has been applied to spherical particles at 
fluid interfaces. We have presented a thorough validation of the method for single particle situations and 
have shown that a transition from a bijel to a Pickering emulsion can be found by varying the particle 
concentration, the particle's contact angle, or the volume ratio of the solvents. Further, we investigated 
the temporal evolution of the droplet/domain growth in emerging Pickering emulsions and bijels. 

Modelling colloidal particles as perfect spheres is a strong simplification of systems appearing in 
nature. There, the particles are generally not spherical, but might show geometrical distortions or fully 
anisotropic shapes, as is, for example, common for clay particles. As a first step to investigate the impact 
of particle anisotropy on the adsorption and stabilization properties, this paper focuses on ellipsoidal 
particles. In addition to the properties of spheres adsorbed at an interface, in the case of anisotropic 
ellipsoidal particles the orientation becomes important and the process of adsorption is in this case 
more complex [11]. Furthermore the anisotropy of the ellipsoids leads in general to a deformation of 
the interface. However, an adsorbed sphere or ellipsoid with a contact angle 6 = 90° does not deform 
the interface in absence of an external potential such as gravitation. For multiple particles capillary 
interactions, which depend on the distance and the orientation of the particles, become relevant 11211 and 
orientational discontinuous phase transitions of the particles can be found |fl3}]. Experimentally it was 
shown that the number of ellipsoidal particles required to stabilize a fluid-fluid interface decreases with 
increasing particle aspect ratio and that a tip-to-tip arrangement is dominant 11411 . 

The remaining sections are organized as follows: In section|2]the simulation method (lattice Boltz- 
mann combined with molecular dynamics) is illustrated. Since studying particle-stabilized emulsions 
demands an exceptional amount of computing resources we focus on specific implementation details of 
our simulation code in section [3] In particular, we highlight specifically code improvements that allow 
to harness the power of massively parallel supercomputers, such as the Blue Gene/P system JUGENE 
at Julich Supercomputing Centre with its ability to run up to 294912 MPI (Message Passing Interface) 
tasks in pai^allel. The following section reports on simulations of single particle adsorption of ellipsoidal 
particles and the formation of bijels and Pickering emulsions. Finally, we conclude in section [51 



2. Simulation method 

The lattice Boltzmann method is a very successful tool for modelling fluids in science and engi- 
neering. Compared to traditional Navier-Stokes solvers, the method allows an easy implementation of 
complex boundary conditions and — due to the high degree of locality of the algorithm — is well suited 
for the implementation onparallel supercomputers. For a thorough introduction to the lattice Boltzmann 
method we refer to Ref. ||6|]. The method is based on a discretized version of the Boltzmann equation 

fi\x + Ci,t + l)=fi{x,t)+n'i{x,t), (1) 
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where f[{x,t) is the single-paiticle distribution function for fluid component c after discretization in 
space X and time t with a discrete set of lattice velocities c, and 



^ _i;-(x.o-r(P'(x,o,u^(v)) 

T 

is the Bhatnagar-Gross-Krook (BGK) collision operator. ff^{p'^,u'^) is the equilibrium distribution func- 
tion and T is the relaxation time. We use a three-dimensional lattice and a D3Q19 implementation (/ = 
1, . . . , 19). From Eq. ([T]), the Navier-Stokes equations can be recovered with density p' (x,?) = L,/)'^(x,f) 
and velocity u' = Y.ifi^il?'^ i^i the low-compressibility and low Mach number limit. If further fluid 
species c' with a single-paiticle distribution function ff {x,t) are to be modeled, the inter-species inter- 
action force 

P'(x,0 = -»F^'(x,0£^.c'I*I'''(x',0(x'-x) , (3) 

c' x' 

with a monotonous weight function ^"^{x,t) for the effective mass is calculated locally according to the 
approach by Shan and Chen and incorporated into the collision term Q-l in Eq. ([T]) |19|]. In our case, the 
coupling strength ^^c' is negative in order to obtain de-mixing and the sum over x' runs over all sites 
separated from x by one of the discrete c;. Colloidal particles are discretized on the lattice and coupled 
to both fluid species by means of a moving bounce-back boundary condition ifisl. if x is part of the 
surface of a colloid then Eq. ([T|) for adjacent fluid sites x + c; is replaced with 

f[{x + Ci,t + l)=ff\x + Ci,t)+n';{x + Ci,t)+C, (4) 

where C = — r-p''(x + c,-,f)usurf • c,- is a linear function of the local particle surface velocity Usurf and 
the direction / is defined via c; = — c^. a^ . and Cg are constants of the D3Q19 lattice. The particle 
configuration is evolved in time solving Newton's equation in the spirit of classical molecular dynamics 
simulations. As the total momentum should be conserved, an additional force 

F{t) = {2ff{x + Ci,t)+C)cj (5) 

acting on the particle is needed to compensate for the momentum change of the fluid caused by Eq. (|4]l. 
The potential between the particles is a Hertz potential which approximates a hard core potential and 
has the following form for two spheres with the same radius R [16J: 

<\>H = KHi2R-rij)'2 for r,y < 2/?. (6) 

rtj is the distance between the two sphere centers and Kh the force constant. For the simulations which 
are discussed later in this text a value of Kh = 100 is used. In the next step the potential is generalized to 
the case of ellipsoids with the parallel radius Rp and the orthogonal radius Ro by following the method 
which was applied by Berne and Pechukas for the case of an intermolecular potential llTI] . We define 
o = 2R and 8 = KhO^ and extend a and 8 to the anisotropic case so that 

8(6,-,6;) = — = and a(6,-,6;,f,';) = — ^ _ =, (7) 

V 2V i+x6,6; l-x6,6; ^ 

with a = 2Ro, X = g2 I d2 . and 6, the orientation vector of particle /. The scaled potential can be written 

as 

^H{oi,dj,rij) = e{6i,bj)^H i ^'J . x ) ■ (8) 
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Figure 1: Relation between the particle color Ap and the contact angle 0. A linear relation is found: = 243.2°Ap + 90°. 
Inset: Definition of the contact angle for an ellipsoidal particle adsorbed at an interface between two fluids. 

is a dimensionless function which takes the specific form of the potential form into account. In 
addition to adding the direct interaction described by the Hertz potential we correct for the limited de- 
scription of hydrodynamics when two particles come very close by means of a lubrication correction. 
If the number of lattice points between two particles is sufficient, the lattice Boltzmann algorithm re- 
produces the correct lubrication force automatically. The error that occurs if the flow is not sufficiently 
resolved can be corrected by 

3nfiR\ f 1 1 \ 

F,7 = ^r,.(r,.(u,. - u,)) - - J (9) 

in the case of two spheres with radius R. We choose a cut-off at = | and u; is the velocity of particle 
/. This equation is generalized to ellipsoids in a similar way as the Hertz potential using Eq. (Q. 

The force in Eq. ([3]) also includes interactions between lattice nodes outside of particles with a lattice 
node inside a particle. To calculate these interactions the lattice nodes x in the outer shell of the particle 
are filled with fluid densities 

P«.w = P' = ^ L P'(x + , t) (10) 

'NP 

corresponding to the average over the Nnp non-particle lattice nodes adjacent to x in directions i^p. We 
consider a system of two immiscible fluids, which we call red and blue fluid, and in which particles are 
suspended. By defining a parameter Ap, the particle color, we are able to control the interaction between 
the particle surface and the two fluids and thus control the contact angle as it is defined in the inset of 
Fig.[T] If Ap has a positive value, we add it to the red fluid component as p^g„ = p' + Ap. Otherwise we 
add its absolute value to the blue fluid as p^^„, = p ^ + |Ap|. In Fig. [T]it is depicted that the dependence 
of the contact angle on the particle color can be fitted by the linear relation 

= 243.2°Ap + 9O°, (11) 

where the slope depends on the actual simulation parameters. For a more detailed description of our 



simulation algorithm the reader is referred to Ref. llOll . 



3. Implementation 

Development of our simulation code LB3D started in 1999 as a parallel LB solver capable of describ- 
ing systems of up to three fluid species coupled by Shan and Chen's aforementioned approach fl fisll . 
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Figure 2: Strong scaling of LB3D on the Blue Gene/P before and after the optimizations. |(a)| relates to a system with only 
one fluid component so that the effect of matching or mismatching topologies of network and domain decomposition can be 
examined better. |(b)| refers to a system with two fluid species and suspended particles as they are of interest in this paper. The 
absolute execution times for small core counts did not change significantly. 



In 2008, a — since then severely extended — parallel molecular dynamics code (l^ was integrated into 
LB3D using a common 3-dimensional spatial decomposition scheme. It is employed here to implement 
the colloidal particles. 

The strong locality of the LB equation Eq. ([T|l as well as the short-range interaction forces rele- 
vant in colloidal systems generally allow for an efficient parallelization. However, while several authors 
presented highly efficient combined single-component LB -suspension codes for the massively paral- 
lel simulation of blood flows on IBM Blue Gene/P systems |[20j, |2l|], it still is not trivial to achieve 
good scalability on this platform. For example, the network employed for MPI point-to-point commu- 
nication provides direct links only between nearest neighbors in a three-dimensional torus. Allowing 
MPI_Cart_create ( ) to reorder process ranks and manually choosing a domain decomposition that fits 
the known hardware topology can therefore increase the performance significantly at high degrees of par- 
allelism. Fig. 2(a)] demonstrates this on the basis of strong scaling speedup for a system of 1024^ x 2048 
lattice sites carrying only one fluid species and no particles. Consequently, only Eq. ([T|l demands sig- 
nificant computational effort, which makes possible communication bottlenecks more visible. On the 
other hand, JUGENE, the IBM Blue Gene/P system at Julich Supercomputing Centre, consisting of 
294912 cores, brings to light serial parts of a code which would stay undetected at lower core counts. 
Fig. |2(b)| visualizes the effect of parallelizing a loop over the global number of colloidal particles present 
in earlier versions of our particle-fluid coupling routines which — in this strong scaling benchmark — had 
a visible effect only when scaling beyond 8192 cores. The test system of 1024^ x 2048 lattice nodes 
contains two fluid species and 4112 895 uniformly distributed particles of spherical shape with a radius 
of 5 lattice units at a volume concentration of 20 % and is thus very close to the Pickering systems of 
interest. When running on 262 144 cores, our two-component LB-suspension solver achieves an average 
of 8.97 X 10^ lattice updates per second (LUPS) in total at a decomposition of 32 x 16^ LB sites per 
core. This corresponds to 3.42 x 10"* LUPS per core compared to 3.66 x 10^ LUPS per core for the 
minimum core count of 2048 or a relative parallel efficiency of 93.5%. On average, for 262 144 cores, 
41.8% of the computing time is spent on the Shan-Chen force evaluation Eq. Q, followed by 29.6% 
for the particle-fluid coupling and 18.7% for the LB equation Eq. ([B. The remaining costs are mainly 
caused by communication. 

While the benchmai^ks above relate to the pure time evolution, actual production runs require check- 
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Figure 3: (a) Snapshots of the adsorption of an elhpsoid with aspect ratio m = 2 and a contact angle = 90° for a starting 
angle of (|) = 0.0057°. Every snapshot is related to a square symbol in Fig. |3(b)| (b) Adsorption diagram for an ellipsoid with a 
parallel diameter Rp = 12, an orthogonal diameter Ro = 6, and a contact angle = 90°. z is the distance between the particle 
center and the interface, (|) is the angle between the ellipsoid main axis and the interface normal. The adsorption trajectories 
are presented by the dashed lines. The solid lines depict where the particle would touch an infinitely thin and undeformed 
interface. The circles correspond to the stable ((|) = ^) and the metastable ((|) = 0) equilibrium point. The squares correspond 
to the snapshots in Fig.|3(a)| 



pointing and the output of physical observables to disk. Using parallel HDF5 output, we succeed in 
storing fluid density fields of 4.6GB size for a system of 1024^ x 1 152 lattice sites and 454508 particles 
in 29s (on average) when employing the whole system (294912 cores). This corresponds to the time 
required to simulate about 100 LB steps and is acceptable since output is required less than once per 100 
to 1000 time steps. 

Finally, it is important to be awai^e of possible peculiarities of the MPI implementations found on 
large supercomputing systems. Due to their relatively small size in memory we collect and write par- 
ticle configuration data serially on the root process. When running on 131072 cores of JUGENE, we 
encounter a speedup of 77 when employing the specially optimized MPI_Allgatherv ( ) compared to 
the more intuitive MPI_Gatherv ( ) , at the cost of requiring one receive buffer per task. 

4. Results and discussion 

After having discussed the implementation of our program code some of the results obtained by this 
program will be presented. We study a system of two immiscible fluids (called red and blue fluid in 
the following) stabilized by prolate ellipsoidal particles. A first step to understand fundamental physical 
properties of such systems is the study of the behavior of a single ellipsoidal particle at a flat fluid-fluid 
interface. If the particle touches the interface the particle is adsorbed to it. The binding energy of the 
particle is usually of the order lO^^^r so that the adsorption is an irreversible process. This binding 
energy arises from the fact that an adsorbed particle reduces the interfacial area between two fluids with 
repulsive interactions. The adsorption was studied in Ref. [ 1 1] using a free energy method. There, the 
free energy is calculated for all particle orientations and particle center-interface distances. Then, the 
equilibrium state is found by minimizing the free energy function. In contrast to our simulations, in 
the approach given in Ref. fllO the interface is infinitely thin and undeformable, while we take a finite 
interface thickness and a interface deformability into account. 

To simulate a single-paiticle adsorption we use a cubic volume with a side length of 64 lattice units 
with two lamellae of different fluids where each lamella has a thickness of 31 lattice units. The boundary 
conditions are periodic in the two directions parallel to the fluid-fluid interface. Impenetrable walls are 
installed in the third direction at z = 1 and z = 64 with a thickness of 1 lattice unit each. To keep 
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Figure 4: Two snapshots of emulsions stabilized by particles are shown for an aspect ratio m = 2, a volume concentration of 
C ~ 0.2 and a contact angle = 90°. |(a)| Bijel for a fluid ratio of 1 : 1. [(b)] Pickering emulsion for a fluid ratio of 5 : 2. 

the interface thin enough the fluid-fluid interaction parameter between red and blue (Eq. ([S])) is set to 
gbr = 0.1. Snapshots of the adsorption for a single ellipsoidal particle with a contact angle & = 90° and 
an aspect ratio m = ^ = 2 are shown in Fig. |3(a)| for a starting angle (|) = 0.0057°. Rp and Rg are the 
ellipsoidal radii parallel and orthogonal to the symetrie axis. It can be observed that at the beginning of 
the simulation (first 2600 time steps in the present example) the particle moves in the direction of the 
interface without changing its orientation significantly. Thereafter (t = 2600... 6500 steps) the orientation 
changes from (|) ?a 0° to (|) w 90° and the particle reaches its equilibrium position. 

Each snapshot in Fig. |3(a)| is related to a square symbol in Fig. |3(b)[ Here, z is the distance between 
the particle center and the interface, (j) is the angle between the ellipsoid main axis and the interface 
normal. The solid lines correspond to the points where the particle just would touch an infinitely thin 
and undeformed interface. The dashed lines correspond to the adsorption trajectories and the two equi- 
libria are illustrated by circles. The upper circle (at (|) = 90°) corresponds to the stable point being in 
relationship with the global free energy minimum whereas the lower one (at (|) = 0°) shows a metastable 
point. The value of the z coordinate of the stable and metastable point depends on the contact angle. In 
the example shown above both points have the value z = 0. If the adsorption trajectory of the particle 
starts with (|)(? = 0) = 0° it reaches the metastable point. For any other starting angle {(^{t = 0) / 0°) 
the adsorption trajectory of the particle ends in the stable point. The adsorption lines approach attrac- 
tor lines. For initial angles much lai^ger than (|) = 0° there are two attractor lines, one on each side. 
The adsorption trajectories coming from the respective sides approach this attractor line. There is a 
third attractor line at z = 0. The adsorption trajectories starting at very small angles approach this at- 
tractor Une. Qualitatively comparable theoretical calculations of the adsorption have been presented in 
Ref. flT]. Quantitative comparisons are not easily possible since the theoretical approach assumes an 
infinitely thin and undeformable interface, while our diffuse interface simulations take the interfacial 
deformability into account. 

In the following we consider emulsions stabilized by a large number of anisotropic colloidal parti- 
cles. To simulate a bulk system we use a cubic volume with a side length of 256 lattice units and periodic 
boundary conditions. To obtain a phase separation the fluid-fluid interaction parameter between red and 
blue (see Eq. ([3])) is set to gbr = 0.08. The initial density distributions of the two fluids ai^e chosen to be 
random. 

As stated in the introduction we distinguish between two different phases for these emulsions, the 
Pickering emulsion and the bicontinuous interfacially jammed emulsion gel (bijel). The bijel (see 
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Figure 5: Left plot: Phase diagram demonstrating the transition from a bijel to a Pickering emulsion. The contact angle and 
the ratio between the two fluids are varied. The squares show the configurations which lead to a bijel whereas the circles denote 
a Pickering emulsion. Right plot: Average domain size L in equilibrium versus the volume concentration of the particles C for 
a bijel with particles of neutral wettability (0 = 90°). The plot compares the data for spheres (m = 1, cubic symbols) fioll with 
ellipsoids (here: m = 2, circles). The behavior can be described by L = ^ + (see solid line (m = 2) and dashed line (m = 1)), 
where a and b are fit parameters. 



Fig. |4(a)| ) consists of two continuous phases whereas the Pickering emulsion |[l], 0] (see Fig. |4(b)[ ) con- 
sists of droplets of one fluid immersed in a second fluid phase. Bijels and Pickering emulsions stabilized 
by spherical particles have been investigated in Ref. jlOll . In this paper simulation results for ellipsoids 
with an aspect ratio of m = 2 are presented. The choice of the control parameters decides into which of 
the two phases the system evolves. The left plot in Fig.[5]shows the transition from a bijel to a Pickering 
emulsion for an aspect ratio of m = 2 and a particle volume concentration of C « 0.2. The two control 
parameters used for the study of the phase transition are the fluid ratio and the contact angle 0. The 
squares show the configurations which lead to a bijel whereas the circles denote a Pickering emulsion. 
If the amount of the two fluids present in the simulation is equal or not too different (e.g. a ratio of 
4 : 3) we find a bijel for all considered contact angles. However, if the fluid ratio is increased we find 
a Pickering emulsion. For intermediate fluid ratios the obtained phase depends on the chosen contact 
angle. For example for a ratio of 9 : 5 we get a bijel for a contact angle of 90° and a Pickering emulsion 
for all higher values of 0. 

To characterize an emulsion the time dependent lateral domain size L{t) = i(Lv(?) +Ly{t) +L^{t)) 
is calculated. Its Cartesian components (/ = x,y,z) are defined as 

MO^^, (*?('» -5^^. (.2) 



{kf{t)) is the second-order moment of the three-dimensional structure function 5(k,f) = ^|<l)k(OI- 4*' = 
$ — ($) is the Fourier transform of the fluctuations of the order parameter $ = p,- — ph- The right plot in 
Fig.|5]shows the dependence of the domain size on the volume concentration of particles C for different 
concentrations between C = 0.08 and C = 0.24 for an emulsion with a fluid ratio of 1 : 1 and = 90°. 

The results for anisotropic particles (here: m = 2) and spheres (m = 1) are compared. As can be 
seen in the left inset of Fig. [5] these parameters for the case of m = 2 lead to a bijel. For too small 
volume concentrations such as C 0.04 there are not sufficiently many particles available to stabilize 
the emulsion. The system finally reaches a fully demixed state. The simulation data is displayed by 



circles for m = 2 and squares for m = 1 (see Uvil ). L decreases with an increasing value of C. An 



increasing number of particles stabilizing the emulsion leads to an increase of the total interfacial area 
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and thus to smaller domain sizes. The results can be fitted by the equation L{C) = ^ + b |[lol] (see 
solid (m = 2) and dashed (m =1) lines in the right plot in Fig. O with values b ^ 1 1.08 for m = 2 and 
b ^ 10.85 for m = 1. The offset b is due to the finite system size in our simulation and is approximately 
independent of the particle shape. The factor a depends on the shape of the particles. For spheres we 
find a value of t? « 3.86 whereas for m = 2 we find a 2.20. The right plot in Fig. |5] shows that the 
domain size L is larger for spherical particles than for the case of anisotropic particles. If a spherical 
particle with a radius R and a contact angle = 90° is adsorbed it reduces the interfacial area by an 
amount of A(m = 1) = 2nR^. If an anisotropic ellipsoid with aspect ratio m and the same volume as the 
sphere is adsorbed and reaches the stable point (see results of the single-particle adsorption (Fig. |3(a)[ )) 
the interfacial area is reduced by an amount of A{m) = InRpRg = Inm^l^R^. For the example of m = 2 
we obtain a reduction of the interface area of A{m = 2) ^ 2nl.26R^. Thus, if the particle anisotropy 
increases, the occupied area A (m) increases. This means that if the particle number is kept constant and 
the anisotropy is increased, the particles can stabilize a larger interface which leads to smaller domain 
sizes. As the structure of a bijel interface is quite complicated and the interfaces are generally curved it 
is not possible to utilize the relations given above to derive an exact dependence of the domain size L on 
the particle aspect ratio m. However, the qualitative difference of using spheres or ellipsoids as emulsion 
stabilizers can be well understood. 

5. Conclusion and outlook 

We have demonstrated simulations of anisotropic ellipsoidal particles stabilizing fluid-fluid inter- 
faces based on a combined multicomponent lattice Boltzmann and molecular dynamics approach. We 
provided a description of the code implementation and its recent improvements enabling our code for 
efficient large-scale simulations of Pickering/bijel systems on hundreds of thousands of cores on an IBM 
Blue Gene/P system while sampling physical observables such as density fields or colloidal particle con- 
figurations with sufficient temporal resolution. After basic studies of adsorption trajectories for single 
ellipsoids we demonstrated that ellipsoidal particles can lead to a transition between bijels and Pickering 
emulsions depending on contact angle, particle volume concentration or fluid ratio. Further, we demon- 
strated that the average lateral domain size of bijels depends on the particle concentration and can be 
fitted by a simple 1/C relation where the pref actor depends on the particle aspect ratio demonstrating 
that particles with m > 1 are more efficient emulsion stabilizers than spheres. 
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